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Abstract 

We study a master equation system modelling a population dynamics problem in 
a lattice. The problem is the calculation of the minimum size of a refuge that can 
protect a population from hostile external conditions, the so called critical patch 
size problem. We analize both cases in which the particles are considered fermions 
and bosons and show using exact analitical methods that, while the Fermi-Dirac 
statistics leads to certain extinction for any refuge size, the Bose-Eistein statistics 
allows survival even for the minimal refuge. 
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1 Introduction 



Biological systems are known to be some of the most complex systems found in 
Nature. Lately, they have been receiving great atention from the mathematical 
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and physical sciences [1], that try to give fundamental explanations to the 
phenomenology found in the field. One of the most important and changelling 
problems in ecology is to understand the mechanisms that lead one population 
to extinction, interesting from both the fundamental and applied point of view. 
Due to the complexity of the problem the most of the approaches have been 
phenomelogical, like the use of stochastic differential equations [2] or reaction- 
diffusion equations [3]. Another consequence of this complexity is the lack 
of an axiomatic theory describing this kind of systems. However, it has been 
noted that there is a deep connection between ecological and reaction-diffusion 
processes [4]; a connection that may be exploted to get a more fundamental 
understanding of ecology. 

The problem under consideration is the so called critical patch size problem, 
already classic in the mathematical literature. It was first considered by Skel- 
1am in his influencing article of 1951 [3], where he studied a population living 
in a finite refuge where the conditions are good for life, but with hostile con- 
ditions outside. He showed that there is a critical size of the refuge such that 
if the actual refuge is smaller than it the population get extincted. He used a 
reaction-diffusion equation approach to that end, that considers a continuum 
population living in a continuum space. But since the population is actually 
discrete and finite, we will use a different approach here, say, a master equa- 
tion approach, that considers a discrete population living in a discrete space. 
This type of models have been studied via Monte Carlo simulations and mean- 
field aproximations [5,6], but here, instead, we will perform exact analitical 
calculations. It is worth noting that Skellam's results have been improved re- 
cently in a model that incorporates the internal fluctuations effects due to 
the discreteness of the population [4] (see also [7] for a related problem on 
extinction). This model, as well as Skellam's theory, deals with a continuum 
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space and suppose every individual of the population punctual. The present 
model deals with a lattice, that may be interpreted as a correction to the in- 
finitesimal size of the individuals, while it represents the lost of the continuum 
space. This kind of analitical calculations are interesting also as a guide for 
numerical simulations that take place, of course, on a lattice. 

In this work we will get some insight into the extinction issue via two different 
models of reacting-diffusing particles: in the first one the particles will obey 
the Fermi-Dirac statistics while in the second they will obey the Bose- Einstein 
one. The biological interpretation is straightforward: the Fermi-Dirac statistics 
impose a maximum in the number of particles (individuals) per lattice site, 
i. e., a strict carrying capacity of the medium. In the case of Bose-Einstein 
statistics there is no such limit, even when some "crowded" realizations of 
the stochastic process could be highly unlikely. It is easy to see that it is 
absolutely necessary to choose either statistics before studying a reaction- 
diffusion problem, be it numerical or theoretical. 

At first look it seems that the statistics chosen has no relevance in this problem 
since the usual approach is to consider the system near extinction [3,4,8]. In- 
deed, the standard procedure is to linearize the corresponding equation around 
the null population value, this is, to analize the low occupation number aprox- 
imation. And since we are studying a system that seems to be controlled by 
its low occupation number state, this suggests that the requirements on the 
occupation number (the particle statistics) will have no effect on the possi- 
ble extinction event. It is actually a question of theoretical interest if a sys- 
tem of reaction-diffusion bosons behave similarly to fermions in such systems 
dominated by a low occupation number state (see, for instance, the discus- 
sion in [9]). We will show that this is not the case. Actually, for the case of 
fermionic particles, an arbitrarily well-adapted population will get extincted 



3 



in finite time due to a rare fiuctuation event. This fact is very interesting 
also from the point of view of the stochastic modeUing, since rare fluctua- 
tions are some of worst understood issues that appear both in theory and 
experiments, and whose understanding and control is very important for the 
pure science and applications respectively [10]. Also, the development of new 
techniques to study reaction-diffusion fermionic flows is a very important sub- 
ject that has proven itself very interesting for both points of view pure and 
applied [11,12,13]. In the case of bosonic particles the situation changes com- 
pletely and the population will be able to survive even for the minimal refuge 
size. This implies that the importance of the statistics chosen for modeling 
the population is complete. 

In order to proof this statements we will use an interesting analogy with 
quantum mechanics, already common in the study of reaction-diffusion prob- 
lems [9,14] and that has proved itself useful when apphed to biological sys- 
tems [4,15]. 

2 The Fermionic Model 

Our master equation will model a population of random walkers living in a 
flnite lattice, that will be able to reproduce, but will not neither compite for 
the nutrients nor die inside the refuge; death is only allowed outside. This 
model can be thought as unreal, but we will study it because it overestimates 
the possibilites of survival: if the population get extincted in this model, it will 
get extincted if we considered death and/or competion. Every site of the lattice 
will be able to keep at most A'' individuals, indicating the finite resources of the 
medium to maintain a population. If we choose a low value of N, like = 1, 
we are modelling a starving population, but because we can choose arbitrarily 
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high values of N, we can model arbitrarily well adapted populations. 

Let us begin with the simplest case: A three sites one-dimensional chain un- 
dergoing the following series of reactions. The first site (top to the left) allows 
diffusion to the right 

A + 0^0 + A (1) 
at rate D, reproduction to the right 

A + $^A + A (2) 

at rate a, and death 

A^$ (3) 

at rate 7. The third site (top to the right) is symmetric to this one and allows 
death, reproduction to the left, and diffusion to the left at the corresponding 
rates. The middle site allows diffusion and reproduction to both the left and 
the right, but it is considered absolutely safe for the population, so no death 
reaction is allowed. Additional reactions of death and competion can be added, 
but they will only rend the extinction more likely. We will further consider a 
two-state chain: every site of the lattice can be either empty or occupied by a 
single individual. This system can be described via a master equation of the 
form: 

dP(t) ^ 

= nmj - ^)m - w{i ^ j)pmi (4) 

where Pi{t) describes the probability of being in the state |i) at time t, and 
W{j ~^ i) is the transition rate from the state |j) to the state It is worth 
pointing out that the algebraic properties of the master equation guarantee 
that its solution consist of a linear combination of terms with a decaying 
exponential time-dependence, and so will always show a stable approach to 
some steady state [16,17]. 

The steady states of this equation are given by the condition ~ 0, and 
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in our particular case this condition reduces to the hnear algebraic equation: 



M-V^O, 



(5) 



where M is the 8x8 matrix 
/ 



-2j a 2a a 



-T D a a 00 
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7 D 



-T 
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where T — a-\- D -\- ^ and V is the vector 



(|111) , |110) , |101) , |011) , |100) , |010) , |001) , |000))* 



(6) 



that must fuUfil the normalization condition 



|111) + |110) + |101) + |011) + |100) + |010) + |001) + |000) = 1. (7) 



In our notation 1 stands for an occupied site and for an empty site, and 
\iik) for the probability of being in such a state. The left hand side of Eq.(5) 
can be thought as a linear map; this way our problem reduces to calculate its 
kernel. It is easy to see, by performing the matricial product and solving the 
corresponding linear system, that, provided a > 0, -D > 0, and 7 > 0, the 
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kernel is generated by the linear span of the vector: 

(0,0,0,0,0,0,0,1)*. (8) 

This is also the only element of the kernel that satifies the normalization 
condition Eq.(7), which means that it is the unique steady state solution of 
the master equation. Physically, this means that, no matter with which initial 
condition are we starting, the system will be in the state |000) in the infinite 
time limit. That is, extinction is certain for long times. 

This result is not surprising in view of the continuum space calculations of 
critical patch sizes, where one can show that a patch shorter than the critical 
one leads to certain extinction, for a critical patch size strictly greater than 
zero [3] . In this discrete version of the problem, we have chosen the minimum 
possible size of the refuge, so this result could be expected a priori. However, 
we would like to see what happens in situations with different number of sites. 
This way we will see if it is possible to define a critical patch number, that is, 
to determine what is the minimum number of sites being part of a refuge that 
rends it effective to prevent an extinction in the infinite time limit. 

Let us consider a one-dimensional two-state finite lattice with L + 2 sites. The 
first site (top to the left) allows diffusion to the right, reproduction to the right 
and death at the rates defined below. The {L + 2) —th site (top to the right) is 
defined as symmetric to the first one. The L central sites allow reproduction 
and diffusion both to the left and the right, but not death. Note that allowing 
only two sites to be out of the refuge rends more difficult the extinction than 
if we consider more, since the individuals outside can only diffuse inside. It is 
interesting to compare this with the continuum case, where it is considered an 
infinite space out of the refuge at both boundaries. 

In this case we can construct a master equation like in the former one, and 
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write the corresponding steady state condition generalizing Eq.(5). This will 
be again a linear map given this time by a 2^^+^^ x 2^^+^^ matrix. Note that 
all the elements in the 2*^^+^'* — th column of this matrix are identically zero, 
reflecting the fact that |0...0) is an absorbing state: once we arrive at this state, 
we will stay there forever. This is the mathematical expression of the physical 
fact that no population can be created from nothing. This means that |0...0) 
always solves the master equation, something that should appear as obvious 
if we realize that it can be chosen as initial condition, and therefore we will 
stay there for every t > 0. For L = 1 we have also shown that the vacuum 
state is the unique steady state solution of the master equation, implying that 
extinction is certain. We can establish therefore that the nonuniqueness of 
steady state solutions of the master equation is a necessary condition to avoid 
extinction. This is equivalent to search for situations in which the dimension 
of the kernel of the linear map has a dimension greater than or equal to two. 
Provided that the determinant of the matrix is always zero (due to the fact that 
the smallest possible kernel is one-dimensional), we are looking for situations 
in which the determinants of the minors of the matrix are all zero. Fix L > 1; 
in this case we have 2^-^+^^ x 2*^^+^^ or less algebraic equations of at most 
[(2(^+2) - 1) X (2(-'^+^) -l)]-th order, which should be all zero in order to get a 
two-dimensional kernel. Suppose that all this equations but one are identically 
zero, this way we arrive at a situation that is not closer to extinction than the 
real one (we will consider the case of all equations identically zero below) . Fix 
a,D > 0; now we have a single algebraic equation in the variable 7 of at most 
[(2(^+2) _ 1) X (2(^+2) - 1)] -th order, that has at most (2^^+^) - 1) x (2(^+2) - 1) 
solutions, as a consequence of the fundamental theorem of algebra. Note that 
7 = is a solution, since in this case we will arrive certainly to the state 
in the infinite time limit, provided that the initial condition is other than the 
vacuum state (in this case we can consider the null solution as "unstable"). 
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Suppose that there are < (2(^+2) - 1) x (2(^+2) - 1) solutions that are 
real and positive, and denote 7* the smallest of these solutions. For 7 = we 
have survival, as we have shown, and let us start now varying continuously 
the value of 7 from zero to upper values. Because the set of solutions of an 
algebraic equation is countable and finite, there is an open interval (0,7*) of 
values of 7 that cause certain extinction. That is, for values of 7 greater than 
zero but infinitely close to it, the system is driven to extinction. As we vary 
continuously 7 we get certain extinction for an infinite number of values till 
we arrive to 7*, when we get again survival. This means that we get survival 
for higher values of the death rate, while extinction is certain for lower values 
of it, something that is absurd. We must necessarily conclude that there is not 
such 7*, or what is the same, extinction is certain for every 7 > 0. Note that in 
the case of more than one equation not identically zero, we can reproduce the 
same argument defining 7* as the smallest common root of all the equations 
different from zero. 

Let us examine now the case in which all the equations are identically zero. 
Suppose that there is a value of L, say L, that fullfils this property. In this 
case, there should be at least another solution different from the vacuum state, 
at which we should arrive if we start from the appropiate initial condition, no 
matter what the values of the parameters are. Let fix 7, D > 0, and a = 0; 
in this case is easy to see that we will get certain extinction in the infinite 
time limit for every initial condition, a contradiction. We must necessarily 
conclude that there is not such Z, or what is the same, for every L there is at 
least one equation that is not identically zero. This means that, provided that 
a,^,D > 0, we get certain extinction in the limit t 00. 

Generalizations of this problem are straightforward. Suppose that we have an 
hypercubic (i-dimensional lattice with L sites per side. Suppose that we allow 
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A'" individuals per site, and diffusion and reproduction to the first neighbours, 
and, in the case oi N > 1, also on-site reproduction. Suppose further that we 
allow the death reaction to take place in just the boundary sites. A similar 
argument as the above one will lead us again to the same conclusion: ex- 
tinction is certain. Even in more complicated cases, like when the occupation 
number depends on the lattice site N = N{i), where this dependence may 
be deterministic or stochastic (quenched disorder), we get the same result via 
an equivalent argument. We can further claim that, because the probability 
of survival is identically zero in the steady state, all the realizations of the 
stochastic process should reach the vacuum state in finite time. 

The physical implications of this fact are very important. We can consider arbi- 
trarily well adapted populations by considering arbitrarily high (but bounded) 
values of N and appropiate initial conditions. We can further choose high val- 
ues of the reproduction rate and low values of the death rate; no matter of 
this the population will finally get extincted. This means that the dynamics 
of this system is dominated by a rare fiuctuation that always appears in finite 
time and kills the whole population. 



3 The Bosonic Model 



In this section we will analize the same model but with bosonic particles for 
contrast. The Malthusian population under consideration will undergo birth, 
A ^ A + A, a.t rate a, at every site, death, A ^ 0, at rate 7, only at the 
boundary sites, and diffusing coupling between first neighbours. To do the 
analysis we will need a very different methodology. The master equation in 
this case reads 
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dP{{n,}- ^) _ ^ ( ^ ^[(^^ + _ 1^ + 1^ . ^) 



+a[(ni - 1)P(.., rii - 1, t) - niP{..., m, t)]) , (9) 

where {e} denotes the set of first neighbours of the site i and i is not at the 
boundary. In the case i is a boundary site the master equation is 



= E f^E[K + l)P{■.■,n^ - l,ne + 1, 

i ^ {e} 

-niP{...,ni,ne, + 
(7[(nj - rij - 1, t) - niP{..., n^, t)] + 

^[(m + 1)P(..., m + 1, t) - niP{..., m, t)]) . (10) 

The analytical treatment of these equations comes as follows. We can map 
this master equation description of the system into a quantum field-theoretic 
problem. This connection was first proposed by Doi [18] and was deeply gen- 
eralized in subsequent works [14]. We can write this theory in terms of the 
second-quantized bosonic operators: 

[4, ttj] = Sij, [ai, ttj] = 0, [al, a]] = 0, |0) = 0, (11) 

whose effect is to create or to destroy particles at the corresponding lattice 
site: 

al\...,ni,...) = \...,ni + l,...) , (12) 



ai \...,ni, ...) = rii |...,nj - 1, ...) , (13) 
where we have defined the states as: 

\{n^}) = m^lr\o)■ (14) 
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Thus we can define the time-dependent state vector as: 



\m) = T.PiM;t)\M), (15) 

{Ui} 

and claim that it obeys the imaginary time Schrodinger equation 

±mt))^-Hm), (16) 

with the hamiltonian 



H^JZi-^T. al(«e - ai) + a[l - al]ala, . (17) 

i \ {e} J 

Note that we recover Eq.(9) if we substitute Eq.(15) and Eq.(17) in Eq.(16). 
In the case of the hnear chain, this second quantized theory is equivalent to 
the stochastic differential equation 

da — D{ai+i + ai- 1 - 2ai)dt + aadt + \/2aadW{t), (18) 

where dW{t) denotes the increments of a Wiener process, and the correct 
interpretation of the equation is Ito. This connection is explained in [14], 
where it is also proved that the first moment of the stochastic process a{t) is 
the first moment of the population. This implies 

^ = D{{ai+^) + (a,_i) - 2{ai)) + (7(a,), (19) 
at 

in any site outside the boundary. If we repeat the calculation we get 

d{ai) 



dt 



^D({ae)-{ai)) + (a-^){ai), (20) 



when i is at the boundary and e denotes the first neighbour of i. This allows 
us to write the exact equation for the mean values of the population in the 
minimal chain (L=3): 

- = M.a, (21) 
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where a — (oi, 02, 03)* and M is the 3x3 matrix 



(7-7-D D 



D a -2D D 



\ 







D 



We choose as initial condition a — (0, 1, 0), because is the minimal initial con- 
dition that keeps one particle under the refuge. The solution can be calculated 
straightforwardly, but since the expresions are cumbersome we will only report 
here the solution at the boundary, at, — ai — a^, 



ab{t) = -D 



exp 



i (-3D - 7 - y/9n' - 2D7 + 7^ + 2a) t 



+D- 



exp 



V9D^ - 2D7 + 7^ 
i (-3D - 7 + V9D^ - 2D7 + 7^ + 2a) t 



(22) 



V9D^ - 2D7 + 7^ 

It is easy to see that the solution is increasing in time at some site if and only 
if it is increasing in time at every site, thus it is enough to study the solution 
at one site, for instance, ah{t). And ab{t) is increasing if and only if 

1 



(J>- (^D + 7 - -y9D2^^2D^ 



r 



(23) 



So we have seen that in the bosonic model the minimal chain allows survival 
for the appropiate parameter values. 



4 Conclusions 



In this work, we have seen that the statistics of the reacting-diffusing particles 
modelling one ecological population has a fundamental importance when we 
study the issue of extinction. It could be thought a priori that this would 



13 



not be the case since whenever we are close to extinction, due to the low 
number of surviving particles, restrictions on the occupation number would 
not be important. We have shown that this is actually not the case, and that 
while fermions are certainly driven to extinction for any refuge size, bosons 
can survive even in the minimal refuge. 

We can also conclude that the risk of extinction can only be defined like a 
probability of survival in a finite interval of time in the case of fermions. This 
definition has already been used in biology [5,19], and it rules out the possibil- 
ity of use those mathematical models that suppose an absorbing boundary for 
enough high values of population. We have shown that this supposition can 
not be assumed for enough long times, and we have put in a more rigorous 
footing the definition first heuristically elucidated by the biologists. 

But, in the case of bosons, the mathematical treatment for finite times is not 
necessary, because the infinite time limit aproximation is enough to predict 
and separate surviving populations to those driven to extinction. Finally, as a 
remark, we would like to underline the methodological character of this work, 
that tries to clarify the mathematical properties of the models commonly used 
to understand ecology, rather than give a realistic picture of the extinction 
problem for one concrete population. 

Very interesting mathematical problems are still to be solved, for instance, to 
find the analytical expression of the time dependent solution of the master 
equation in the case of the one-dimensional fermionic chain with L + 2 sites. 
This type of problem is usually projected onto a spin chain problem expressed 
via Pauli operators, that turns out to be exactly integrable in some cases [20]. 
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